Asymmetric spreading in highly advective, disordered environments. 
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Spreading of bacteria in a highly advective, disordered environment is examined. Predictions 
of super-diffusive spreading for a simplified reaction-diffusion equation are tested. Concentration 
profiles display anomalous growth and super-diffusive spreading. A perturbation analysis yields 
a crossover time between diffusive and super-diffusive behavior. The time's dependence on the 
convection velocity and disorder is tested. Like the simplified equation, the full linear reaction- 
diffusion equation displays super-diffusive spreading perpendicular to the convection, ffowever, 
for mean positive growth rates the full nonlinear reaction-diffusion equation produces symmetric 
spreading with a Fisher wavefront, whereas net negative growth rates cause an asymmetry, with a 
slower wavefront velocity perpendicular to the convection. 

PACS numbers: 87.23.Cc,87.18.-h,05.40.-a 
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The study of population growth is an integral part of 
the biological sciences. Recently the behavior of micro- 
bial species, such as bacteria, has enjoyed much mathe- 
matical analysis due to the formation of intricate equi- 
librium patterns P, Q ■ Due to environmental and health 
concerns the behaviors in driven systems is also of inter- 
est. In this regime the dynamics are typically described 
using reaction-diffusion equations which may involve a 
number of species and their interactions . While much 
work has been done on these types of systems, the inclu- 
sion of disorder in the environment has received limited 
attention. This disorder may manifest itself in a variety 
of ways, from spatial variations in the available food or 
in the presence of poisons to random diffusion constants. 

In particular, a reaction diffusion equation with spa- 
tially varying growth factors may take the following form, 

dtc(x, t) = DV^c{x, t)-v- Vc(f , t) 

+ [a + U{x)]c{x,t) -bc^{x,t), (1) 

where the U{x) are spatially random growth rates and 
the convection velocity v drives the population through 
the environment 0]. If both U = and v — then 
Eq. n reduces to the Fisher equation, where the growth 
rate a provides exponential growth which is cut off by 
the nonlinear interaction term b at the system's carrying 
capacity, a/6 0. The linear regime of Eq. ^ around the 
fixed points c = and c = a/b, has been studied in detail 
in Ref . (3(| . Only some limited numerical simulation have 
been performed for the nonlinear case |3| ■ In the limit of 
large convection velocity u, fixed in the y direction, and 
with 6 = 0, a substitution of the form 
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allows one to obtain a simplified form of Eq. ^ 

vdtW{x, t) = D\l'^W{x, t) + U{x, t)W{x, t). (3) 

Here y is relabeled as t and the directions perpendicu- 
lar to y as Q . As the substitution (Eq. |2) contains 
the exponential growth and diffusion in the y direction, 
the function W{x, t) describes the cross section perpen- 
dicular to the convection for a population at the 'time' 
t — y/v. Interestingly, this simplified equation has the 
form of an imaginary time Schrodinger equation with a 
random, fluctuating potential. Additionally it is directly 
connected to the problem of directed polymers in random 
media [HI . It has been shown that Eq. 13 predicts super- 
diffusive growth for the long time, large distance behav- 
ior, with an exponent of 2/3 in one spatial dimension |^. 
This exponent has been reproduced numerically by ex- 
amining the averaged mean squared displacement of the 
optimal path (lowest energy path) in directed polymers 
[3 ■ Through a detailed examination of the full concentra- 
tion profiles of Eq. |S1 this report examines the behavior 
of this simplified equation in the context of population 
growth with the goal of obtaining a better understanding 
of the full equation's (Eq. ^ behavior in both the linear 
(6 ~ 0) and nonlinear (6 > 0) regimes. First a perturba- 
tion analysis of Eq.|2|yields a crossover time dividing pure 
diffusion and the super diffusive behavior. Afterwards, 
one dimensional numerical simulations describe concen- 
tration profiles and anomalous growth of Eq. 13 and test 
the predictions for the diffusion exponent and crossover 
time. 

To obtain a perturbation expansion for Eq.|3J first note 
that it is an initial value problem. Thus, in the spirit of 
Ref. 0, one applies a Fourier-Laplace transform. 
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FIG. 1: Diagrammatic representation of the Fourier-Laplace 
transformed, simplified equation (Eq. |SJ|. Part (a) gives the 
exact equation while part (b) displays the disorder averaged 
perturbation series to one loop order. 



Equation 13 then takes the form, 
Wik,uj) = vGo{k,io)W{k,0) 



xU{q,na)W{q,ni,), 



(5) 



with the abbreviations J d'^q = ^^^^^ , j dt = dt, 
and f dVl,, = where y = a or 6 where 

J y Jy— zoo Itti ^ 

Go{k,uj) — {vuj + Dk^)^^ is the free propagator and W 
denotes taking only the Fourier transform. In obtain- 
ing this form, the Bromwich integral giving the inverse 
Laplace transform was used. 

A graphical representation of Eq. and its second or- 
der, disorder averaged expansion are shown in Fig. Q^a) 
and (b) respectively. As U has zero mean, upon tak- 
ing the disorder average the first order term drops out. 
Hence determining to one loop order the renormalized 
propagator G_r, defined as W{k,uj) = vGR{k,uj)W{k,0), 
requires calculating only the second order term. For 
a uniform distribution of width A the correlator is 
t)U{x',t') = ^l%5'^{x - x')5{t - 1') where ly is the 
lattice constant for the y direction. With this fact and 
some straightforward contour integrations, the renormal- 
ized propagator of Fig. ^^b) becomes. 



Gfl(fc,cj) = Gn{k,u) + 



Gi{k,oj). 
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where Sd is the surface area of a d-dimcnsional sphere of 
unit radius. Expanding Go and G^ for A; — > yields 



24dv^uj 
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When the second term on the right side of Eq. [7| be- 
comes on the order of one, then pure diffusion is no longer 
the dominant term. The crossover time is proportional to 
the value of u;^^ at this point. From Eq. |7|the crossover 
time T is given by, 



T : 
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The crossover time depends on both the velocity of the 
flow as well as the width of the distribution of random 
growth rates. When the velocity increases, the system 
is pushed through the random environment before it has 
time to experience the fluctuations, seeing an effectively 
averaged environment. Hence it makes sense that T is 
increased by higher velocities. On the other hand, in- 
creasing the width of the random distribution of growth 
rates creates optimal paths in the system which have 
larger effective growth rates. This causes pure diffusion 
to break down earlier, hence the inverse dependence with 
the crossover time. 

Numerical simulations of Eq. |3were performed in one 
spatial dimension using a Runge-Kutta technique 0- A 
Gaussian initial condition of unit variance was centered 
on a lattice of 20000 sites. This was large enough to 
insure that the boundaries were never encountered by 
the concentration. The random growth rates depend on 
time, and so must be updated during the simulation. To 
provide equal time and space lattice constants, a Runge- 
Kutta step size of 0.1 was chosen and the growth rates 
were updated after every ten time steps. Lastly, the con- 
centrations were normalized after every time step. 

The initial concentration profile mimics the inocula- 
tion of a medium with an initial bacterial sample. With 
no disorder present the Gaussian shape would persist, 
with a variance increasing in time. However, the dis- 
order destroys this by providing particularly favorable 
paths along which growth may occur. In Fig. El concen- 
tration profiles are shown for two different times after 
inoculation. The top two plots contain profiles for a sin- 
gle disorder realization. For the short time there is very 
little deviation from the purely diffusive situation. On 
the other hand, the large time concentration proflle devi- 
ates greatly, with large spikes developing. These regions 
correspond to the end of a path in the {x, t) space which 
had particularly favorable growth rates and thus resulted 
in a much larger population of bacteria than would be ex- 
pected from a homogeneous environment. For the long 
times, these paths may end with increasing distance from 
the original starting position effectively shifting the mean 
position of the concentration from the starting position. 
There may also be several competing paths that have 
nearly the same effective growth rate resulting in sev- 
eral concentration spikes in the proflle. These behaviors 
do not occur in the homogeneous case as the spreading 
Gaussian profile always remains centered on the starting 
point and symmetric about that point. Upon perform- 
ing a disorder average, as seen in the bottom of Fig. El 
the large, off-center peaks result in average concentra- 
tion profiles whose widths increase faster then the purely 
diffusive case. 

The concentration profiles in Fig. [3 were normalized to 
have a clear comparison with the case of pure diffusion. 
As no growth terms are present in Eq. (31 when U = 0, 
the latter case remains normalized. This breaks down in 
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FIG. 2: Concentration profiles for one dimensional spread- 
ing according to Eq. |3 normalized by the total population. 
The top plots are for a single disorder realization while the 
bottom plots are averaged over 1000 random configurations. 
A Gaussian of unit variance was used as the initial profile. 
The solid lines are profiles for a disordered environment with 
V = 4.0 and A = 1.0 and the dashed lines show the purely dif- 
fusive case. A maximal and several characteristic error bars 
are shown on the averaged distributions. 



the former case however. Even though [/ has a zero aver- 
age, fluctuations in the growth factors lead to anomalous 
growth. Although the effective growth is relatively small, 
corresponding to a = 0.00223/^t in Eq. ^for the system 
shown in Fig. |2 at the longest times (t = SOOOOOZj) the 
effects are profound as the total population becomes on 
the order of e^^"*'. Clearly a diverging bacterial density is 
unphysical. It emerges here because the nonlinear death 
term has been dropped in the simplified equation, Eq.|31 
A simple argument for the appearance of this anoma- 
lous growth lies in the asymmetry between the growth 
and death processes. In particular, consider a small con- 
centration of bacteria present in a favorable environment, 
[/ > 0. This concentration will grow exponentially in 
time and will spread additional concentration to neigh- 
boring areas via diffusion. On the other hand, in an 
unfavorable environment for growth, with C/ < 0, the 
local population will experience an exponential decay in 
concentration. While this decreases the total local pop- 
ulation, the decrease imposed in neighboring areas is not 
the opposite of the growth case. The asymmetry lies in 
the fact that the local (and total) population cannot fall 
below zero. One cannot have a negative number of or- 
ganisms. Hence the difference in concentration between 
neighboring sites, proportional to the rate of transport, 
is smaller with C/ < leading to a reduced rate of pop- 
ulation loss compared to the gain in population when 
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FIG. 3: Super-diffusive behaviors for the simplified equation 
in one dimension: (a) Concentration width vs. time averaged 
over 1000 disorder configurations with v = 4.0 and A = 4.0 
(solid line) and the purely diffusive case (dashed line). The 
two curves begin with a power-law of 0.5, but for long times 
the disordered case switches to a power-law of 0.647 ± 0.001, 
shown as a dot-dash line, (b) Crossover times for variation 
of A with V — 4.0 (squares) and variation of v with A — 1.0 
(circles). The solid line shows a power law of — 1.97±0.05 and 
the dashed line a power law of 2.03 ± 0.10 as fit to the data 
varying A and v respectively. Each point is averaged over 
3000 random configurations with statistical error bars on the 
order of the symbol size. 



U > 0. Therefore the asymmetry of the diffusion pro- 
cess in the random environment leads to an effectively 
increased growth rate. 

To measure the diffusion exponent, the width of the 
disorder averaged concentration profile, as shown in Fig. 
12 was measured as a function of time. The resulting 
curve is shown in Fig. Efa) . For pure diffusion the profile 
width grows as as one expects. The disordered case 
is different, with two regions of clearly different power- 
law behaviors. For small times the disordered width fol- 
lows the behavior of the purely diffusive case. However, 
as the time becomes large, the disordered case deviates 
from pure diffusion and instead grows with a power-law 
exponent equal to 0.647 ± 0.001. This super-diffusive 
behavior is in good agreement with the exponent value 
2/3 that has been previously predicted Q. As explained 
above, this super-diffusive behavior arises due to the ap- 
pearance of optimal growth rates that deviate far from 
the center of the population. Upon averaging, these shift 
concentration from the center of the profile resulting in 
a width wider then the diffusive case. 

A clear crossover to super-diffusive behavior is seen in 
Fig. Eta). However, the location of this crossover depends 
on the simulation parameters. For the simulations, the 
predicted crossover time, Eq.|Hl becomes T = 247rt>^/A^. 
Compared to the visually apparent crossover point in Fig. 
|3{a), the prediction of T « 75 is roughly an order of mag- 
nitude too small. This should not be completely unex- 
pected as Eq.lHlreally describes the time where departure 
from purely diffusive behavior begins. At this point the 
width is growing super diffusively, but does not saturate 
at the full exponent until i 3> T. 
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FIG. 4: Two-dimensional asymmetric spreading for the full 
reaction-diffusion equation (Eq. with convection along the 
vertical axis: (a) disorder averaged concentration maps for 
the linear case, b — 0, show super-diffusive spreading perpen- 
dicular to the convection velocity for A = 2, whereas (b), the 
nonlinear case 6 > with a < and A = 3, has a wavefront 
that propagates more slowly in the direction perpendicular 
to the convection, resulting in the opposite asymmetry. For 
comparison, the dashed lines show the homogeneous, linear 
case in (a), and a circle, such as is obtained for the nonlinear 
case with fe > and a > 0, in (b). 

The crossover time's dependence on the velocity and 
random width of the growth rates is shown in Fig.UJb). 
These times were obtained at the point where the differ- 
ence in width between the disordered and purely diffu- 
sive case was equal to that of a baseline case {v — 4 and 
A = 1) at roughly the crossover time predicted by Eq. 
|S1 t = 1200. The crossover times agree very well with 
the quadratic behavior, T ~ (w/A)^, predicted by Eq. |S1 
From Fig.l^Jb) the variation with A at fixed v behaves as 
T ^ /\-i-97±o.05 g^j-^^ ^Yie variation of v at fixed A results 
in a power-law T ^ ^2.03±o.io^ 

As the simplified equation describes the cross section 
of Eq. ^ perpendicular to the convection, it implies that 
Eq. ^ with 6 = should exhibit super-diffusive behavior 
in that direction. Indeed, as seen in Fig.^Ja), the concen- 
tration contours for a two-dimensional simulation of this 
linear case shows contours of width equal to the homo- 
geneous case in the direction parallel to the convection, 
but spreading faster perpendicular to it. 

As noted above, the linear case is not physical for long 
times due to unrealistic organism densities. The nonlin- 
ear case with 5 > presents a much different spreading 
picture. For a > a symmetric Fisher wave develops. 
This symmetry should not be completely unexpected. In 
the linear case, the mechanism for the enhanced spread- 
ing perpendicular to the convection was the disorder av- 
eraging of the large, asymmetric concentration spikes, 
such as shown in Fig. |21 Since these spikes are cut off 
by the carrying capacity, this behavior is absent in the 
nonlinear case and spreading is symmetric. The spread- 
ing is still enhanced, however, as the wavefront velocity 
increases with increasing disorder. In particular, out- 
side the wavefront the linear regime applies and enhanced 
growth is found. Thus, one may argue ,S] that the growth 



rate in the Fisher velocity expression should be replaced 
by the real growth rate, giving the wavefront velocity 



2^/{a,s{A)+a)D. (9) 

Here acs is the effective growth rate of the corresponding 
linear problem which depends on the disorder strength A. 
Numerical simulations of Eq. ^ for a range of disorders 
find excellent agreement with this wavefront velocity 

The wavefront velocity expression, Eq. has an im- 
portant implication. Namely, attempting to poison or 
destroy a colony of organisms, by applying a < 0, may 
fail if the disorder creates sufficiently enhanced growth, 
Ucs > —a. Even more interesting, the resulting wave- 
front is asymmetric but in the opposite manner to the 
above linear case. Figure 0fb) shows the wavefront ob- 
tained from a numerical simulation of Eq. ^ with & > 
but a < 0. The direction parallel to the convection has a 
wavefront velocity that follows Eq.|51 but the perpendicu- 
lar wavefront velocity is smaller, resulting in an asymmet- 
ric droplet. Qualitatively, the smaller wavefront velocity 
arises due to the loss of optimal growth paths. The net 
negative growth rate does not allow paths passing regions 
of random, negative growth rates which were previously 
possible due to the additional positive growth factor. The 
parallel direction is unaffected, because detours around 
these lost paths may take place on each side of the lost 
paths. However, if the lost path occurs on the edge of 
the growing droplet, only paths nearer the droplet center 
remain, resulting in a reduced spreading speed. 
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